


Institutional Archive of the Naval Postgraduate School 


Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1996-03 


Parametrics of submarine dynamic stability in 
the vertical plane 


Papanikolaou, Stavros I. 


Monterey, California. Naval Postgraduate School 


http://ndl.handle.net/10945/32197 


Downloaded from NPS Archive: Calhoun 


| Calhoun is the Naval Postgraduate School's public access digital repository for 
D U DLEY research materials and institutional publications created by the NPS community. 
get Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS'‘s first 
KNOX appointed — and published — scholarly author. 





LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 


http://www.nps.edu/library Monterey, California USA 93943 





NAVAL POSTGRADUATE SCHOOL 
Monterey, California 





THESIS 


PARAMETRICS OF SUBMARINE DYNAMIC 










STABILITY IN THE VERTICAL PLANE 
by 
Stavros I. Papanikolaou 


March, 1996 





Thesis Advisor: 


Approved for public release; distribution is unlimited. 


Fotis A. Papoulias 

















| REPORT DOCUMENTATION PAGE Form Approved OMB No. 0704 | 


Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instruction, searching existing data 






sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other 
aspect of this collection of information, including suggestions for reducing this burden, to Washington headquarters Services, Directorate for Information Operations and 

Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188) 
ashington DC 20503. 















2. REPORT DATE 
March 1996 


3. REPORT TYPE AND DATES COVERED 
Master’s Thesis 


1. AGENCY USE ONLY (Leave blank) 












. TITLE AND SUBTITLE Parametrics of Submarine Dynamic Stability {| 5. FUNDING NUMBERS 
in the Vertical Plane. 
6. AUTHOR(S): Stavros I. Papanikolaou 


. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 
Naval Postgraduate School 
Monterey CA 93943-5000 





8. PERFORMING ORGANIZATION 
REPORT NUMBER 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


. SPONSORING/MONITORING AGENCY NAME(S) AND 
ADDRESS(ES) 









11. SUPPLEMENTARY NOTES The views expressed in this thesis are those of the author and do not reflect the official policy 
or position of the Department of Defense or the U.S. Government. 









12a. DISTRIBUTION/AVAILABILITY STATEMENT Approved for 12b. DISTRIBUTION CODE 


public release; distribution unlimited 





13. ABSTRACT (maximum 200 words) The problem of dynamic stability of submersible vehicles in the dive plane is examined 
utilizing bifurcation techniques. The primary mechanism of loss of stability is identified in the form of generic Hopf 
bifurcations to periodic solutions. Stability of the resulting limit cycles is established using center manifold approximations 
and integral averaging. The hydrodynamic coefficients are calculated using existing semi-empirical methods. Parametric 
studies are performed with varying vehicle geometric properties. The methods described in this work could suggest ways to 
enlarge the submerged operational envelope of a vehicle early in the design phase. 

















14. SUBJECT TERMS 15. NUMBER OF PAGES 
84 


16. PRICE CODE 











Submarine stability, Bifurcations, Periodic solutions 





18. SECURITY 19. SECURITY 20. LIMITATION OF ABSTRACT 
CLASSIFICATION CLASSIFICATION UL 
OF THIS PAGE OF ABSTRACT 





17. SECURITY 
CLASSIFICATION OF 
REPORT 





Unclassified Unclassified 


mer ht RS a+ nn Py SSS, SSS SSS SSS SP chy arya ls ticishatiter % 


NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std. 239-18 


Unclassified 


il 














Approved for public release; distribution is unlimited. 


PARAMETRICS OF SUBMARINE DYNAMIC STABILITY IN THE VERTICAL 
PLANE 


Stavros I. Papanikolaou 
Lieutenant Jounior Grade, Hellenic Navy 
B.S., Hellenic Naval Academy, 1989 


Submitted in partial fulfillment 
of the requirements for the degree of 


MASTER OF SCIENCE IN MECHANICAL ENGINEERING 


from the 


NAVAL POSTGRADUATE SCHOOL 
March 1996 


Author: 
. Papanikolaou 





Approved by: 
- Papoulias, Thesis Advisor 


Terry/R/ Mcnelley, Chairm 
Department of Mechanical Engineering 


ili 


1V 











ABSTRACT 


The problem of dynamic stability of submersible vehicles in the dive plane is examined 
utilizing bifurcation techniques. The primary mechanism of loss of stability is identified in the 
form of generic Hopf bifurcations to periodic solutions. Stability of the resulting limit cycles is 
established using center manifold approximations and integral averaging. The hydrodynamic 
coefficients are calculated using existing semi-empirical methods. Parametric studies are 
performed with varying vehicle geometric properties. The methods described in this work could 
suggest ways to enlarge the submerged operational envelope of a vehicle early in the design 


phase. 
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I. INTRODUCTION 


A. PROBLEM OVERVIEW 


The increasing demands of using submersible vehicles for more complex and 
demanding missions, force us to use a variety of methods, mathematical mod- 
els, and assumptions for the study of their dynamic interactions and responses. 
This study is important in order to enhance vehicle operations. Typically, lin- 
earization of the equations of motion around nominal straight line level flight 
paths along with eigenvalue analysis can be employed (Arentzen and Mandel, 
1960), (Clayton and Bishop, 1982), (Feldman,1987). A simple but efficient 
stability criterion Gv > Q can be obtained where the stability index G, is 
function of the hydrodynamic coefficients in heave and pitch. Values for the 
stability index can be computed by, 


My(Zq +m) 


Gy=1- 
ZuM, 


(1) 


This index is analogous to the familiar stability coefficient for horizontal plane 
maneuvering and can be thought of as a high speed approximation where the 
effect of the metacentric restoring moment is minimal (Papadimitriou,1994). 
If the value of G, is greater than zero, the vehicle is dynamically stable. As it 
has been established in previous studies (Papoulias and Papadimitriou, 1995) | 
though, this is only a sufficient, and rather conservative condition for stability. 
Nevertheless, it is widely used and its value is indicative of vertical plane sta- 


bility for any new design. We should keep in mind, however, that the condition 








G,, < 0 indicates a divergent loss of stability which is quite uncommon in the 
vertical plane. Most modern submarines exhibit a flutter—like instability at 
high speed, which can not be analyzed using the above simplified index. Di- 
vergent motions may develop in combined six degrees of freedom (Papoulias et 
al, 1993) and their occurrence can not be analyzed by a single stability index. 
Previous work (Papadimitriou, 1994) was limited to a single body with fixed 
hydrodynamic coefficients. In this work, we expand by allowing the geometry 


of the body and thus its hydrodynamic properties to vary. 


B. THESIS OUTLINE 


Previous work (Papoulias and Papadimitriou, 1995) analyzed the problem 
of stability of motion with controls fixed in the vertical plane, with partic- 
ular emphasis on the mechanism of loss of stability of straight line motion. 
The closed loop control problem was analyzed in (Papoulias et al, 1995). The 
surge equation was decoupled from heave/pitch through a perturbation series 
approach (Bender and Orszag, 1978). As was established in (Papadimitriou, 
1994) loss of stability occurs in the form of generic bifurcations to periodic 
solutions (Guckenheimer and Holmes, 1983). Taylor expansions and center 
manifold approximations were employed in order to isolate the main nonlinear 
terms that influence system response after the initial loss of stability (Hassard 
and Wan, 1978). Integral averaging was performed in order to combine the 


nonlinear terms into a design stability coefficient (Chow and Mallet—Paret, 














1977). Some difficulties associated with the nonsmoothness of the absolute 
value nonlinearities was dealt with by employing the concept of generalized 
gradient (Clarke, 1983). This was employed as an alternative to the lin- 
ear/cubic approximation typically used in ship roll motion studies (Dalzell, 
1978). The same methodology is applied in this work in order to analyze the 


sensitivity of the results with respect to geometric characteristics of the body. 


Vehicle modeling in this work follows standard notation (Gertler and Ha- 
gen, 1976), (Smith et al, 1978), and numerical results are presented for a family 
of bodies of revolution similar to the DARPA SUBOFF model (Roddy, 1990) 
for which a set of hydrodynamic coefficients and geometric properties is avail- 
able. This parametric study is conducted utilizing existing semi~empirical 
methods for the calculation of hydrodynamic coefficients. The methods are 
based on (Fidler and Smith, 1978), (Humphreys and Watkinson, 1978), (Pe- 
terson, 1980) and have been verified in (Wolkerstorfer, 1995). The effects of 
varying the nose, base, and tail fractions of the body as well its nondimen- 
sional volume to length ratio on the hydrodynamic derivatives were studied in 
(Holmes, 1995) where prediction equations were derived based on curve fitting 
of the results. These hydrodynamic prediction equations are normalized by 
taking the SUBOFF model as a baseline. This model has been experimentally 
validated for angles of attack on the hull between +15 deg., while the constant 
coefficient approximation introduces very little error in time domain simula- 
tions (Tinker, 1978). Unless otherwise mentioned, all results in this work are 


presented in standard dimensionless form with respect to the vehicle length 








L = 4.26 m, and nominal forward speed U = 2.44 m/sec (Papadimitriou, 


1994). 











Il. PROBLEM FORMULATION 


A. EQUATIONS OF MOTION 


In order to obtain the mathematical model the following assumptions, re- 


strictions, and definitions have to be made: 


1. The submersible vehicle motion is restricted in the vertical plane, thus 


the model consists of coupled nonlinear heave and pitch equations. 
2. The coordinate frame is fixed at the vehicle’s geometrical center. 
3. Vehicle is port/starboard symmetric and neutrally buoyant. 


4. Use Newton’s equations of motion in dimensionless form. 


The nonlinear heave and pitch equations become: 


m(w — ug — zag’ — 264) = ZG + ZwwW + Zqq + Zww 


nose 
-Co [ ' b(x)(w — xq)|w — xq| dz, (2) 
Iq + mzg(ut+ wa) — mz¢e(w — ug) = Mag t+ Myw + Mqqt+ Mow 
nose 
+Cp | b(x)(w — xq)|w — xq|x dx 


tai 


—zgpW cosé — zgpW sin 8, (3) 


where tgp = XG — XB, 2GB = 2G — Zp, and the rest of the symbols are based 
on standard notation as shown in Table 1. Without loss of generality we can 


assume that zp = zp = 0, so that xggp = rg and zegp = zg. The cross flow 





integral terms in these equations become very important for high angles of 
attack maneuvering, where they provide the primary motion damping. The 
drag coefficient, Cp, is assumed to be constant throughout the vehicle length 
for simplicity. This does not affect the qualitative properties of the results 


that follow. The vehicle pitch rate is, 
b= ¢. (4) 


Dynamic coupling between surge and heave/pitch is present due to coordinate 
coupling as a result of the nonzero metacentric height. However, it has been 
shown (Papoulias and Papadimitriou, 1995) that this coupling is of higher 


order and does not change the linear and nonlinear results that follow. 


B. HYDRODYNAMIC COEFFICIENTS 


Systematic studies based on semi~empirical methods have resulted in the 
evaluation of hydrodynamic coefficients for a generic body of revolution in 
terms of basic geometric properties. Curve fitting revealed that adequate ac- 


curacy for initial design can be obtained by equations of the form 
Ho = AjF? + AoFpFm+ A3F2 + AaFa 
V 
+AsFy, + Ag+ Az (= — c) ; 
where Hc denotes a given coefficient in its standard nondimensional form, V 


the underwater volume of the body, L its nominal length, F,, the nose fraction, 


and F,,, the mid—body fraction. The regression coefficients A; are presented 
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(xp, 2p) | body fixed coordinates of vehicle center of buoyancy 
2G, zq) | body fixed coordinates of vehicle center of gravity 


Table 1: Nomenclature 
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ia, | 0.0003 | +0.0040 | +0.0027 | -0.0012 | —0.0045 | +0.0008 | —0.1590. 
it, 
My [ =0.0031 


Table 2: Regression coefficients A; 
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Ag 






in Table 2. Z; was assumed constant since the semi-empirical techniques 
failed to compute a reliable value. Basic geometric definitions for the body 
are presented in Figure 1. The constant C’ is approximately 8 x 10-3 and is 
the nominal value for the volumetric coefficient. These expressions are for a 
body of revolution without appendages and assume parabolic nose, parallel 
mid—body, and conical tail (Holmes, 1995). Typical ranges of applicability for 
these regression formulas are 0.05 to 0.25 for F,, 0.40 to 0.60 for F,, and 6.0 to 
10.0 for V/L°. Sample results for the above hydrodynamic coefficients versus 


the nose and mid—body fraction ratios are presented in Figures 2 through 8. 


C. DEGREE OF STABILITY 


The degree of stability is defined as the largest real part of all eigenvalues of 
the linearized system of equations (2), (3), and (4). Positive values indicate an 
unstable system while negative values show stability of forward motion. The 
degree of stability versus rgg for constant forward speed u = 0.0 and different 


values of zgp is shown in Figures 9 through 12. Based on these results we can 
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Figure 1: Geometric definitions 





Figure 2: Hydrodynamic coefficient Mj versus Fy, and Fm 





Figure 3: Hydrodynamic coefficient M4 versus Fp, and Fm 
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Figure 4: Hydrodynamic coefficient Z,, versus F, and Fm 
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Figure 5: Hydrodynamic coefficient M, versus F,, and Fm 
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Figure 6: Hydrodynamic coefficient Z, versus F, and Fm 
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Figure 8: Hydrodynamic coefficient Z versus F, and F'm 
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Degree of stability 





-0.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01 
Xg 


Figure 9: Degree of stability for u = 0.5, varying zqgp, Fn = 0.3, and Fy, = 0.6 
draw the following conclusions: 


1. In all cases the vehicle is dynamically more stable as the metacentric 


height zgg is increased. 


2. In all cases the vehicle is dynamically less stable as the separation between 


the centers of gravity and buoyancy is reduced in absolute value. 


3. For constant F,, increasing values of F,, result in less stable vehicles. 


This means that a longer tail is beneficial for stability of motion, as 


expected. 


4. The same conclusion holds for constant mid—body ratio F,, and varying 


nose ratios F,,. 
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Figure 11: Degree of stability for u = 0.5, varying zg¢p, Fn = 0.3, and Fm = 0.4 
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Degree of stability 
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Figure 12: Degree of stability for u = 0.5, varying zgp, Fn = 0.1, and Fy, = 0.6 


Corresponding results for constant z¢p = 0.015 and varying forward speeds 
u are shown in Figures 13 through 16. Similar conclusions as those discussed 


previously hold in these cases with the following exceptions: 


1. For very low forward speeds, the case xg = 0 may be best for stability. 


2. For very low speeds, smaller tails may result in more stable configura- 


tions. 


Combined results for variations in both zgg and u are shown by the mesh 
plots of Figures 17 through 20. The value of z¢g was held constant at 0.015 
for all plots. These figures confirm our previous conclusions by presenting the 


results in more detail. 


Figure 21 shows the degree of stability versus F,, and F,,. Both values 
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Figure 13: Degree of stability for z¢p = 0.015, varying u, F, = 0.3, and Fm = 0.6 
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Figure 14: Degree of stability for zgg = 0.015, varying u, F, = 0.1, and Fm = 0.4 
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Degree of stability 
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Figure 15: Degree of stability for zgg = 0.015, varying u, Fy, = 0.3, and Fm = 0.4 
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Figure 16: Degree of stability for zgg = 0.015, varying u, F, = 0.1, and Fm = 0.6 
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Figure 17: Degree of stability for F, = 0.3 and Fm = 0.6 
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Figure 18: Degree of stability for F, = 0.1 and Fy, = 0.4 
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Figure 19: Degree of stability for F, = 0.3 and Fi, = 0.4 
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Figure 20: Degree of stability for F, = 0.1 and Fy, = 0.6 


1 


of zg and zg were kept constant and equal to 0 and 0.015 respectively. The 
three surfaces shown correspond to values u = 0.4,0.5,0.6. The upper one 
corresponds to u = 0.6 while the lower one to u = 0.4. It can be seen that 
the degree of stability becomes more negative for decreasing u, and, generally 


speaking, for decreasing F,, and Fyn. 


Figure 22 shows the degree of stability versus F, and Fi. Both values of 
forward speed u and zg were kept constant and equal to 0.5 and 0.015 respec- 
tively. The three surfaces shown correspond to values rg = —0.01, 0, +0.01. 
The upper one corresponds to zg = 0.0 while the lower one to zg = +0.01. It 
can be seen that the degree of stability becomes more negative for increasing 


zg in absolute value, and, generally speaking, for decreasing F,, and Fy. 


Figure 23 shows the degree of stability versus F,, and F,,. Both values of for- 
ward speed u and zg were kept constant and equal to 0.5 and 0.0 respectively. 
The three surfaces shown correspond to values zg = +0.005, +0.015, +0.025. 
The upper one corresponds to zg = +0.005 while the lower one to zg = 
+0.0025. It can be seen that the degree of stability becomes more negative for 


increasing zg, and, generally speaking, for decreasing F, and Fim. 


D. CRITICAL SPEED 


The parameter value where the real part of the dominant complex conjugate 
pair of eigenvalues crosses zero defines the point where linear stability is lost. 


This critical point can be computed by considering the characteristic equation 
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Figure 23: Degree of stability versus Fy and F,, for u = 0.5, zg = 0, and zg 


0.005, 0.015, 0.025 


of the system (Papadimitriou, 1994). Routh’s criterion applied to this can be 


solved for the dimensionless weight, 


(5) 


? 


0 


BoC? 
AgDa1 — BoCo1 


3 


W= 


where, 


= Zw) (ZGB COS A —ZGB sin 60) : 


m 


( 


Co, = 


Zw(ZEB sin Oo — ZGB COS Ao) : 


Do 


It should be mentioned that the effect of the forward speed u is embedded into 


the definition for the dimensionless vehicle weight W through, 


(6) 





The value of the critical speed u, can then be evaluated from (5) and (6). 
Typical results are presented in Figures 24 through 27. A family of critical 
speeds, u,, is shown versus rg with zg as the parameter of the curves. These 
results were obtained for a nose fraction F,, = 0.1,0.3 and mid—body fraction 
Fim = 0.4,0.6. The volumetric coefficient was kept at nominal for all results. 
Vertical plane motions are stable for forward speeds less than the critical speed. 
It can be seen that stability is increasing with increasing zg while rg = 0 is the 
most conservative condition for stability. Therefore, a vehicle which is stable 
when properly trimmed will remain stable for off-trim conditions. The fact 
that a vehicle with a longer aft-body ought to be dynamically more stable is 
confirmed by comparing the results of Figures 24 and 26 to the results shown 
in Figures 25 and 27 respectively. It can be seen that the corresponding critical 
speeds become smaller, thereby reducing the dynamic stability margin, as the 
nose and mid—body fractions are raised. This trend is consistent for all values 


of xg and zg examined. 


Combined plots of the critical speed versus both zg and zg are shown in 
Figures 28 and 29. Figure 28 presents the surfaces for F, = 0.3 and Fy, = 
0.4,0.5,0.6. The uppper surface corresponds to F,, = 0.4. Figure 29 presents 
the surfaces for F,, = 0.5 and F, = 0.1,0.2,0.3. The upper surface corresponds 


to f= 0.1. 


Combined plots of the critical speed versus both F, and F,, are shown 
in Figures 30 through 32. Figure 30 presents the surface when zg = 0.0125 


and zg = 0. Figure 31 gives us a comparative view keeping z¢ = 0.0125 and 
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Figure 24: Critical speed versus xg for F, = 0.1 and Fy, = 0.4 and different values 
of zg 
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Figure 25: Critical speed versus xg for F, = 0.1 and Fy = 0.6 and different values 
of zg 
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Figure 26: Critical speed versus zg for F, = 0.3 and F,, = 0.4 and different values 
of zc 
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Figure 27: Critical speed versus rq for F, = 0.3 and F,, = 0.6 and different values 
of z¢ 
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Figure 29: Critical speed versus rg and zg for Fy = 0.5 and Fy, = 0.1,0.2,0.3 
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using zg = —0.01,0,+0.01 to plot the surfaces as shown. The lower surface 
corresponds to rg = 0. It can be seen that nonzero rg increases the range 
of stability, while the general trend is to increase stability as both F,, and 
Fr, become smaller. A similar plot for zg = 0 and for three values of zg, 
0.005, 0.010, and 0.025 is shown in Figure 32. The lower surface corresponds 
to zg = 0.005 and the higher one to z = 0.025. It can be seen that the 
metacentric height has by far the greatest effect on dynamic stability, while 


the effects of hull geometry are smaller. 


For comparison, a plot of the classical stability coefficient G, from equation 
(1), is shown in Figure 33. The different curves correspond to various mid- 
body fractions, while the nose fraction is kept constant. It can be seen that G, 
is negative throughout. Therefore, it would have predicted an unstable vehicle 
for all ranges of the parameters, which is of course incorrect. Furthermore, G, 
becomes less negative as F,, is increased, which would suggest that dynamic 
stability is increased as the aft—body length is decreased. This is also a false 
conclusion. As we pointed out in the introduction, the classical stability index 


G, should be used with extreme caution. 
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Ili. BIFURCATION ANALYSIS 





A. INTRODUCTION 


The nonlinear bifurcation analysis is based on the general methodology 
used in (Papadimitriou, 1994). The fundamental equations are reproduced 
here for completeness of the presentation. The nonlinear heave/pitch equations 


of motion (2), (3), and (4) are written in the form, 
6 = q, (7) 
W = aw + 419g + a13(zgBcosé + zggsin A) 


+dy(w, q) ate Cl (w, q) ? (8) 


agiw + aoeq + 023(ZGB cos @ + zqgsin 0) 


nQ. 
I 


+dq(w, q) + co(w,q) ? (9) 


where the various coefficients are functions of the hydrodynamic derivatives 


and mass properties, and J,,, I, are the cross flow integrals. 


The system of equations (7) through (9) is written in the compact form 
x = Ax+g(x), (10) 


where 
x = [0,w, 49] , (11) 
is the three state variables vector, and A is the linearized sytem matrix eval- 


uated at the nominal point x9. The term g(x) contains all nonlinear terms 
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of the equations. Hopf bifurcation analysis can be performed by isolating the 
primary nonlinear terms in g(x). Keeping terms up to third order, we can 


write 


g(x) = g(x) + g(x). (12) 
Using equations (7) through (11), the various terms in (12) can be written as, 
(2) _ 9 


2 
gs ) = (Iy — Mj)mzeq? — (mag + Z4)mzqw9 
+d\)(w,4) , (13) 
gs = —(m— Zy)mzgwqt+(mrgt M,y)mzgq° 


+d) (w, 4) , 
and 
(3) _ 


gy = dQ(w,q) + 


2413(2GB sin Ty — Z2GB COS 99)0° . (14) 


—, 
— 
| 


(3) 
dy’ (w,q) + 
7093(2GB sin aN — ZGB COS 90)0° : 


Expansion in Taylor series of d,,, dy requires expansion of the cross flow inte- 


grals Iy, Iz, which require the Taylor series of 


F(E) = lel. (15) 


This expression can be converted into an analytic function using Dalzell’s 
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approximation (Dalzell, (1978), 


35 €9 


5 
E|é| ¥ T6008 + 48 €, (16) 


which is derived by a least squares fit of an odd series over some assumed range 
of €, namely —& < € < &. This approximation has been extensively used in 
ship roll motion studies and is very useful for its intended purpose. However, 
in the present problem it suffers from the several drawbacks (Papadimitriou, 
1994). Instead of Dalzell’s approximation, we employ the concept of general- 
ized gradient (Clarke, (1983), which is used in the study of control systems 
involving discontinuous or non-smooth functions. In this way we approximate 
the gradient of a non-smooth function at a discontinuity by a map equal to the 
convex closure of the limiting gradients near the discontinuity. In our problem 


we write, 


F(E) = El€ol + 2[Eol(€ — 0) + 
sign(&o)(€ — €0)? + f() , (17) 


as the Taylor series epansion of f(€) near €). The sign function in (17) can be 
approximated by, 
sign(f)) = lim tanh (| (18) 
The quantity y is a small regularization parameter and is used for proper 
normalization of the results. Using (18), we can approximate f(€) in the 
vicinity of &) = 0 by, 
ele oe (19) 
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Since 


frew— 29g, (20) 


we can express the non-smooth cross flow integral terms by, 


C 
LS Z, (Eow® — 3Ew7¢ + 3 Ewa" _ E3q°) 
Y 


C 
= SEs — 3Eqw*¢ + 3E3wq" — E4q°) 


where 


B= [  c'6(2)de, (21) 
tai 


are the moments of the vehicle “waterplane” area. 


Using the previous second and third order Taylor series expansions, equa- 


tion (10) is written in the form, 
x= Ax + g(x) +g) (x) . (22) 


If T is the matrix of eigenvectors of A evaluated at the critical point u = uc, 


the linear change of coordinates, 
x=Tz, z2=T x, (23) 
transforms system (22) into its normal coordinate form, 
a= TOATz4+ To g@(Tz) + T1g9 (Tz) . (24) 


At the Hopf bifurcation point, matrix T~’AT takes the form, 


0 —wo QO 
TAT =| apy 0 ~O:| . 
0 QO p 
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where wo is the imaginary part of the critical pair of eigenvalues, and the 
remaining eigenvalue p is negative. For values of u close to the bifurcation 


poit u,, matrix T~'AT becomes, 


ae —(wo + w’e) 0 
T CAT = | (wo +e) a’e 0 
0 0 pt+p'e 


where ¢« denotes the criticality difference 
E=U—Ue, (25) 
and 


a = derivative of the real part of the critical 
eigenvalue with respect to e , 

w = derivative of the imaginary part of the 

critical eigenvalue with respect to e , 


p = derivative of p with respect toe. 


Due to continuity, the eigevalue p + p’e remains negative for small nonzero 
values of e. Therefore, the coordinate zz corresponds to a negative eigenvalue 
and is asymptotically stable. Center manifold theory predicts that the rela- 
tionship between the critical coordinates z,, z2 and the stable coordinate z3 is 


at least of quadratic order. We can then write zs as, 
23 = 1127 + ayQ2129 + A925 , (26) 
where the coefficients, a;;, in the quadratic center manifold expansion (26) 
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need to be determined. By differentiating equation (26) we obtain, 


23 = 20412121 + 039(2122 ig z129) + 20992922 . (27) 
We substitute 21 = —wozq and Zo = woz, and we obtain 

: 2 2 

23 = 112W02, + 2(a22 = Q31)Woz1 29 — &1QW92%p5 . (28) 


The third equation of (24) is written as, 


23 = pz3t Tg?) (rz)| (29) 


(3,3) 


where terms up to second order have been kept. If we denote the elements of 


T and T7! by, 
T= [my], T=[nig], (30) 
then 
dy 
T1g2)(Tz) = | do | , 
d3 


where expressions for dj, dz, d3, and the coefficients @;; are given in Papadim- 


itriou (1994). 
Equation (29) then becomes 
23 = P23 + dg , (31) 
and substituting (26) and the expression for d3 into (31) we get, 
3 = (pany + ngolo5 + ngalss)z] 
+(por + n3elog + 23336 )2122 
+(porgg + ngolor + ng3é37)29 - (32) 
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Comparing coefficients of (28) and (32) we get a system of linear equations 


which yields the coefficients in the center manifold expansion (26). 


Using the previous Taylor expansions and center manifold approximations, 
we can write the reduced two-dimensional system that describes the center 
manifold flow of (24) in the form, 

Zz) = alezy — (wo + w'e)zq + Fi (21, 22) , 
zg = (wo tw'e)zy + a €z + Fo(z1, z2) , 
where F, Fy are cubic polynomials in z; and zg. 


If we introduce polar coordinates in the form, 
z3= Reos¢d, z=Rsind, 


we can produce an equation describing the rate of change of the radial coor- 
dinate R, 


R=a'eR+ P(6)R°+Q(o)R?. 


This equation contains one variable, R, which is slowly varying in time, and 
another variable, ¢, which is a fast variable. Therefore, it can be averaged 
over one complete cycle in ¢ to produce an equation with constant coefficients 


and similar stability properties, 
R=a'eR+KR°+LR’, 


where 





= 3(8ri1 +113 + raz + 3raq) , 


1 Qn 
L = =| Q(¢)de=0. 


Qn 


Therefore, the averaged equation becomes 


R=a'eR+KR’. (33) 


Equation (33) admits two steady state solutions, one at R = 0 which 


corresponds to the trivial equilibrium solution at zero, and one at 


a 
Ro —_ ~ KS , (34) 


This equilibrium solution corresponds to a periodic solution or limit cycle in 
the cartesian coordinates z 1, zo. For this limit cycle to exist, the quantity Ro 
must be a real number. In our case a’ is always positive, since the system loses 
its stability; i.e., the real part of the critical pair of eigenvalues changes from 
negative to positive, for increasing u. Therefore, existence of these periodic 


solutions depends on the value of kK. Specifically, 


e if K <0, periodic solutions exist for e > 0 or u > ue, and 


e if K > 0, periodic solutions exist for « < 0 or u < Ue. 


The characteristic root of (33) in the vicinity of (34) is 
B= —2a’'e, (35) 


and we can see that 
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Figure 34: Nonlinear stability coefficient versus xg for F, = 0.1, Fm = 0.4, and 
different values of Cp 


e if periodic solutions exist for u > u, they are stable, and 


e if periodic solutions exist for u < u, they are unstable. 


B. RESULTS AND DISCUSSION 


Typical results of the nonlinear stability coefficient K are shown in Figures 
34 through 37. Figure 35 presents a plot of K -y versus zg for zg = 0.015, 
F, = 0.1, Fm = 0.6, and for different values of the quadratic drag coefficient 
Cp. It should be emphasized that the use of K -y is more meaningful than the 
use of K, since it properly accounts for the use of the regularization parameter 


y. Numerical evidence suggests that all curves K -y versus rg converge for 
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Figure 35: Nonlinear stability coefficient versus xg for F;, = 0.1, Fm = 0.6, different 
values of C'p 
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Figure 36: Nonlinear stability coefficient versus zg for F, = 0.3, Fm = 0.4, different 
values of Cp 
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Figure 37: Nonlinear stability coefficient versus xg for F, = 0.3, Fy, = 0.6, different 
values of Cp 


y — 0. For practical purposes, values of 7 smaller than 0.001 produce identi- 
cal results. The results of Figure 8 demonstrate the profound effect that the 
quadratic drag coefficient C'p has on stability of limit cycles. All Hopf bifur- 
cations are supercritical (K < 0), and they become stronger supercritical as 
Cp is increased. It is worth noting that results for Cp = 0 produce subcritical 
behavior, kK > Q, which is clearly incorrect. Thus, neglecting the effects of Cp 
would have produce entirely wrong results in the present problem. Additional 
results show that the bifurcations become stronger supercritical as initial sta- 
bility zg is increased. Figure 34 presents similar results with the only difference 
being the value of mid—body fraction F,, = 0.4. It can be seen that smaller 
Fy, for the same F,, which results in longer body tail, may be beneficial for 


stability in the linear sense but it also generates less supercritical bifurcations. 
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Figure 38: Nonlinear stability coefficient versus F,, for rg = 0, Cp = 0.5, and 
different values of Fy, 


This can probably be attributed to the increased responsiveness of the vehi- 
cle. Figures 36 and 37 show the same results for nose fraction F, = 0.3. It 
should be emphasized, however, that altering the fore and aft body lengths 
might infuence the values of Cp which, as we pointed out, is the single most 


important parameter for the nonlinear nature of the bifurcations. 


Figure 38 shows the nonlinear stability coefficient versus F,, for different 
values of F,, while xg = 0 and Cp = 0.5. It can be seen that smaller F, 
for the same F,,, which results in longer body tail, generates less supercritical 


bifurcations. 


Figure 39 shows tne nonlinear stability coefficient versus F, for different 


values of Fy, while zg = 0 and Cp = 0.5. Again it is clear that longer body 
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Figure 39: Nonlinear stability coefficient versus F, for xg = 0, Cp = 0.5, and 
different values of Fin, 


tail generates less supercritical bifurcations. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


This work presented a comprehensive nonlinear study of straight line sta- 
bility of motion of submersibles in the dive plane under open loop conditions. 
A systematic perturbation analysis demonstrated that the effects of surge in 
heave/pitch are small and can be neglected. Primary loss of stability was 
shown to occur in the form of Hopf bifurcations to periodic solutions. The 
critical speed were instability occurs was computed in terms of metacentric 
height, longitudinal separation of the centers of buoyancy and gravity, and 
the dive plane angle. Analysis of the periodic solutions that resulted from the 
Hopf bifurcations was accomplished through Taylor expansions, up to third 
order, of the equations of motion. A consistent approximation, utilizing the 
generalized gradient, was used to study the non—analytic quadratic cross flow 


integral drag terms. The main results of this study are summarized below: 


1. The critical speed of loss of stability is a monotonically increasing func- 
tion of both vertical and longitudinal LCG/LCB separation. This means 
that a vehicle which is stable when properly trimmed will remain stable 


for off—trim conditions. 


2. Loss of stability occurs always in the form of supercritical Hopf bifurca- 
tions with the generation of stable limit cycles. It was found that this is 


mainly due to the stabilizing effects of the quadratic drag forces. 
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3. Even though the quadratic drag forces do not influence the initial loss of 
stability, they have a significant effect on post—loss of stability stabiliza- 


tion. 


4. In general, longer aft body sections seemed to increase the range of linear 
stability but influence adversely the resulting limit cycles upon the initial 


loss of stability. 


It should be emphasized that the occurrence of supercritical Hopf bifurcations 
is an attribute of the open loop system only. Under closed loop control, it is 
possible to experience either supercritical or strongly subcritical Hopf bifurca- 
tions, as shown in [Papoulias et al (1995)]. The latter are particularly severe 
in practice since self-sustained vehicle oscillations may be initiated prior to 
loss of stability, depending on the level of external excitation or the initial 


conditions. 
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APPENDIX 


The following is a list and description of the computer programs used in this 
thesis. The programs are written in FORTRAN or MATLAB. Complete print- 


outs of the programs follow after the list. 


e CRIT_.0O.M 


MATLAB program for calculating the critical speed for 6 = 0. 


e DSTAB.M 


MATLAB program for calculating the degree of stability. 


e HOPF_O.FOR 
FORTRAN program for evaluation of hopf bifurcation formulas using the 


suboff submarine model. 
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7% Program crit_0.m 
7% Evaluation of critical speed for delta=0 


clear 

rho = 1.94; 

g = 32.2; 

L = 13.9792; 


ndi = 0.5*rho*L”2; 
ndad2 = 0.5*rho*L73; 
nd3 = 0.5*rho*L"4; 
nd4 = 0.5*rho*L"5; 


m = 1556.2363/(g*nd2) ; 
md = 1556.2363/g; 

V = md/rho; 

ZE = 0.005; 


while zg<0.026, 


flagil = QO; 

for Fn = 0.1:0.01:0.35; 
flagi = flagi+l; 
flag2 = <0; 


for Fm = 0.3:0.01:0.6; 


flag2 = flag2+1; 

Fb = 1-Fn-Fn; 

d = ((12*V) ./(pi*L* (3*Fm+2*Fn+Fb))) .70.5; 

re = d/2; 

Vn = (2/3*pi*r.”2*L.*Fn) ; 

Mn = Vn*rho; 

Vm = (pi*r.7~2.*Fm*L) ; 

Mm = Vm*rho; 

Vb = (1/3*pi*r.~2*L.*Fb) ; 

Mb = Vb*rho; 

In = Mn.*(1/5*(r.72+(L*Fn) .72)-(3*L*Fn/8) .72); 

Im = Mm/12.*(3*r.72+(L*Fm) .72); 

Ib = Mb.*(3/5*(r.72/4+(L*Fb) .~2)-(L*Fb/4) .~2) ; 

xcb = pixd.72.*(2*L*Fn. *(L*Fm/2+3*L*Fn/8) ... 
-L*Fb. *(L*Fb/4+L*Fm/2))/(12*V) ; 

Leb = L*(Fn+Fm/2)-xcb; 

Tyd = Int+Im+Ib+ (Mn. *(Lcb-5*L*Fn/8).72)... 


+(Mm.* (Lcob-L*Fm/2-L*Fn).72)... 
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+(Mb.* (Leb-L* (Fn+Fm+Fb/4) ) .*2) ; 
Tym = Iyd/nd4; 


% inputs A1,A2,A3,A4,A5,A6,A7,A8 for each coefficient 

A1=[-0.0641, 0.0277, -0.0314, -0.0003, 0.0002, -0.0002, -0.0031]; 
A2=[-0.1149, 0.0499, -0.0559, 0.0040, 0.0007, -0.0007, -0.0046]; 
A3=[-0.0632, 0.0266, -0.0292, 0.0027, 0.0007, -0.0007, -0.0021]; 
A4=[ 0.0670,-0.0283, 0.0310, -0.0012, -0.0008, .0008, 0.0031]; 
A5=[ 0.0732,-0.0301, 0.0316, -0.0045, -0.0016, .0016, 0.0024]; 
A6=[-0.0263,-0.0056, -0.0091, 0.0006, -0.0144, .0144, -0.0013]; 
A7=[-0.5769,-1.6357, -0.0880, -0.1590, -1.8067, .8067, -0.0808] ; 


Ke Oo O O&O 


% Hydrodynamic coefficient prediction equation 
C1 = 8.023e-3; 
for: i-=.1:7.; 
HCm(€i) = A1Ci)*Fn.°2+A2(i)*Fn.*Fm... 
+A3(i)*Fm.~“2+A4(i)*Fn... 
+A5 (i) *Fmt+A6 (i) +A7 (i) *(V/L73-C1) ; 


zqdot = -6.33e-4; 

HCm(8) = zqdot; 

ratio = (0.5686,-1.4357,-0.2658,0.2675,1.1781,-30.5114,0.8149,1.0]; 
HC = HCm./ratio; 


zqdot = -6.33e-4; 
zwdot = HC(5); 

zq = HC(3); 
Zw = HC(1); 


mqdot = HC(7); 
mwdot = HC(6); 


mq = HC(4); 

mw = HC(2); 

Iratio = 0.92943; 

ly = Iym/Iratio; 
cd = 0.015; 

zb = O/L; 

xudot = —-0.05*m; 

xb = O/L; 

xg = Q; 
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Gv = 1 - mw.*(zqtm) ./(2w.*(mq-m. *xg)) ; 
xgb = xg-xb; 
zgb = 2g-Zb; 
for j = 1:length(zg) 
for i = 1:length(xg) 


theta(i,j) = atan(-xgb(i)./zgb(j)); 


a0 
bO 


cO 
cl 


di 


w 
u0 


(m-zwdot) * (Iy-mqdot) - (nwdot+m*xg (i) )*(zqdot+m*xg(i)) ; 

(-zwdot*m-m*mw-zq*m) *xg(i)+(-m*mqtzwdot*mg-zqdot*mw... 
-zq*mwdot-m*mwdot-ly*zwtmqdot*zw) ; 

—m*zw*xg (i) +mg*zw-zq*mw-m*mw 5 

(-m*xg (i) +zwdot*xg (i) +m*xb-zwdot*xb)*sin(theta(i,j))... 

+(-m*zb-zwdot#zg (j) +zwdot*zb+m*zg (j))*cos(theta(i,j)); 

(zw*xg(i)-zw*xb)*sin(theta(i,j))... 
+(zw*zb-zw*zg(j))*cos(theta(i,j)); 


i,j) = bO*cO/(a0*d1-b0*c1) ; 
(i,j) = (1556.2363/ (nd1*w(i,j)))*.5; 


ucr(flag2,flag1) = u0(i,j); 


end 


end 


end 


end 


Fn 
Fm 


0.1:0.01:0.35; 
0.3:0.01:0.6; 


mesh(Fn,Fm,ucr/8) ,grid 
xlabel (’ Fn’) 
ylabel(’Fm’) 
zlabel(’ucr’) 

hold on 


zg=z_et0.01; 


end 
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% Program dstab.m 


/% Matlab program for calculation the degree of stability 


clear 


clear global 


rho = 1.94; 
gZ = 32.2; 
L = 13.9792; 


ndi = 0.5*rho*L"2; 
nd2 = 0.5*rho*L73; 
nd3 = 0.5*rho*L74; 
nd4 = 0.5*rho*L75; 


m = 1556.2363/(g*nd2) ; 
md = 1556.2363/g; 

V = md/rho; 

fig = 1 


for Fn = 0.1:0.2:0.3, 
for Fm = 0.4:0.2:0.6, 
Fb = 1-Fn-Fn; 
d €(12*V) ./ (pi*L* (3*Fm+2*Fn+Fb))).70.5; 
= d/2; 
Vn = (2/3*pi*r.*2*L.*Fn) ; 
Mn = Vn*rho; 
Vm = (pi*r.72.*Fm*L) ; 
Mm = Vm*rho; 
Vb = (1/3*pi*r.~2*L.*Fb) ; 
Mb = Vb*rho; 
In = Mn.*(1/5*(r.°2+(L*Fn) .72)-(3*L*Fn/8) .*2) ; 
Im = Mm/12.*(3*r.~2+(L*Fm) .72); 
Tb = Mb.*(3/5*(r.~2/4+(L*Fb) .~2)-(L*Fb/4) .*2); 
xcb = pixd.°2.*(2*L*Fn.* (L*Fm/2+3*L*Fn/8)... 
-L*Fb.* (L*Fb/4+L*Fm/2) )/(12*V) ; 
Leb = L*(Fn+Fm/2)-xcb; 
Iyd = Int+Im+Ib+(Mn.*(Lcb-5*L*Fn/8).72)... 
+(Mm.*(Lcb-L*Fm/2-L*Fn).*2)... 
+(Mb.*(Lcb-L* (Fn+Fmt+Fb/4)) .*2); 


ol 








Iym = Iyd/nd4; 


% inputs A1,A2,A3,A4,A5,A6,A7,A8 for each coefficient 

Ai=[-0.0641, 0.0277, -0.0314, -0.0003, 0.0002, —-0.0002, -0.0031]; 
A2=[-0.1149, 0.0499, -0.0559, 0.0040, 0.0007, -0.0007, -0.0046] ; 
A3=[-0.0632, 0.0266, -0.0292, 0.0027, 0.0007, -0.0007, -0.0021]; 
A4=[ 0.0670,-0.0283, 0.0310, -0.0012, -0.0008, 0.0008, 0.0031]; 
A5=[ 0.0732,-0.0301, 0.0316, -0.0045, -0.0016, 0.0016, 0.0024] ; 
A6=[-0.0263,-0.0056, -0.0091, 0.0006, -0.0144, 0.0144, -0.0013] ; 
A7=[-0.5769,-1.6357, -0.0880, -0.1590, -1.8067, 1.8067, -0.0808] ; 


7, Hydrodynamic coefficient prediction equation 
C1 = 8.023e-3; 
for i=1:7, 
HCm(i)=A1 (i) *Fn. 7 2+A2(i) *Fn.*FmtA3(i)*Fm.*2+A4(i)*Fn... 
+A5 (i) *Fm+A6 (i) +A7 (i) * (V/L73-C1) ; 


end 


zqdot = -6.33e-4; 

HCm(8) = zqdot; 

ratio = [0.5686,-1.4357,-0.2658,0.2675,1.1781,-30.5114,0.8149,1.0]; 
HC=HCm./ratio; 


zqdot = -6.33e-4; 
zwdot = HC(5); 
Zq = HC(3); 
zw = HC(1); 


mqdot = HC(7); 
mwdot = HC(6); 


mq = HC(4); 

mw = HC(2); 
Iratio = 0.92943; 

iy = Jym/Iratio; 
ca = 0.015; 


52 





evalsi = eig(A,B); 
degstabi (i, j) 


end 
end 


0.015; 
O/L; 





1556 .2363/ (g*nd2) ; 
-Q.05*m; 


O/L; 


linspace(-0.01,0.01,41); 
8*linspace (0.2, .6,41) ; 


1 - mw.*(zqtm)./(zw.*(mq-m. *xg)) ; 


1556 .2363./(nd1i.*uo.72); 


W; 
xg-xb; 
ZEZ-Zb; 


atan(-xgb./zgb) ; 
1:length(uo) 
1:length (xg) 


[-2*cd 0 0 Q; 
0 Zw (zq+m) 0; 
0 mw (mq-m*xg(i)) (xgb(i)*sin(theta(i))... 

-zgb*cos(theta(i)))*b(j); 

0 0 1 0]; 

[m-xudot 0 m*zg 0; 

0 m-zwdot —(m*xg(i)+zqdot) 0; 

m*zg -(mwdot+m*xg(i)) iy-mqdot 0; 

0 0 0 1); 


figure (fig) 


% no surge coupling 


= max(real(evalsi)); 


mesh(uo/8,xg,degstabi) ,grid 
xlabel(’uo’) 


O3 


ylabel(’xg’) 

zlabel(’degree of stability’) 
fig=figti; 

end 


end 


o4 
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PROGRAM HOPF .FOR 
EVALUATION OF HOPF BIFURCATION FORMULAS 
USING THE SUBOFF SUBMARINE MODEL 


Cd=0.5, ZG=0.015 ,Fn=0.24, Fm=0.52 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION L,IY,MASS,MQDOT,MWDOT,ND1, 
MQ ,MW,K1,K2, 
BETA,GAMA, 
EO ,E1,E2,E3,E4, 
DW1 ,DW2 ,DW3 , DW4, 
DQ1 ,DQ2,DQ3 ,DQ4, 
MASSM,MASSN,MASSB,IB,IM,IN, 
RHO ,CD,RADI , VOLM, VOLN, VOLB,LCB,XCB 


NN DO FP WH DY 


DOUBLE PRECISION M11,M12,M13,M21,M22,M23, 
M31 ,M32,M33, 
N11,N12,N13,N21,N22,N23, 
N31 ,N32,N33, 
L21,L22,L23,L24,L31,L32,L33,1L34, 
L25 ,L26,L27 ,L35,L36 ,L37, 
L21A,L22A,L23A,L24A,L314, 
L32A ,L33A, L34A 

DOUBLE PRECISION LN,LM,LB,FM,FN,FB,KK 


ND FP WD DN & 


DIMENSION A(3,3) ,T(3,3) , TINV(3,3) ,FV1(3) ,IV1(3) , YYY(3, 3) 
DIMENSION WR(3) ,WI(3) , TSAVE(3,3) , TLUD(3,3) , IVLUD(3) 
DIMENSION ASAVE(3,3),A2(3,3) ,XL(55) , BR(55) 

DIMENSION VECO(55) , VEC1(55) , VEC2(55) , VEC3(55) , VEC4(55) 
DIMENSION HCA1(7) ,HCA2(7) ,HCA3(7) ,HCA4(7) ,HCA5 (7) 
DIMENSION HCA6(7) ,HCA7(7) ,HC(8) ,RATIO(7) , SVLUD (3) 


OPEN (20,FILE=’DATA_O.DAT’ ,STATUS=’ NEW’ ) 


00 


mana aq 


& 


WEIGHT= 


RHO = 
DO 8886 
CD = 
G = 
XB = 
DO 8887 
ZG = 
ZB = 
MASS = 
BOY = 
VOLUME= 
DO 8888 
DO 8889 


WRITE 
‘WRITE 
WRITE 
WRITE 


FB 
LN 
LM 
LB 
DIAM = 


WRITE(*, 


RADI 
VOLN 
MASSN 
VOLM 


MASSM = 


VOLB 
MASSB 
IN 





1556 .2363 
13.9792 
1.94 
CDi = 0.40,0.60,0.10 
0.5*CD1*RHO 
32.2 
0.0 
KK = 0.0050,0.0250,0.0050 
KK*L 
0.0 
WEIGHT/G 
WEIGHT 
MASS/RHO 
FN=0.10,0.32,0.10 
FM=0.40,0.62,0.10 


(20,*) °CD =’,CD 
(20,*) ’ZG =’ ,KK 
(20,*) ’FN =’,FN 
(20,*) ?FM =’, FM 


1.0-FN-FM 

L*FN 
L*¥FM 
= L*FB 
SQRT((12.*VOLUME) 
/(3.14159*L* (3. *FM+2 .*FN+FB) )) 


4001) DIAM 
= DIAM/2. 
= (2./3.*3.14159*RADI**2.*L*FN) 
VOLN*RHO 
(3.14159*RADI**2.+*FM*L) 
VOLM*RHO 
= (1./3.*3.14159*RADI**2.*L*FB) 


= VOLB*RHO 


= MASSN*(1./5.* (RADI **2+ (L*FN) **2) 
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LCB 





— (3*L*FN/8) **2) 


= MASSM/12.*(3.*RADI**2.+(L*FM) **2) 
= MASSB*(3./5.* (RADI**2/4. 


+(L*FB) **2)—-(L*FB/4. ) **2) 


= (3.14159*DIAM**2* (2. *L*FN 
* (L*¥FM/2.+3.*L*FN/8.) 
-L*FB* (L*¥FB/4.+L*FM/2.)))/(12.*VOLUME) 
WRITE(*,4001) XCB 


WRITE(*,4001) LCB 

= IN+IM+IB+MASSN* (LCB~5*L*FN/8) **2 
+MASSM* (LCB-L*FM/2-L#*FN) **2 

IY = IY+MASSB* (LCB-L* (FN+FM+FB/4) ) **2 


LY 


= L*(FN+FM/2.)-XCB 


Inputs A1,A2,A3,A4,A5,A6,A7,A8 for each coefficient 


HCA1(1) 
HCA1 (2) 
HCA1 (3) 
HCA1 (4) 
HCA1(5) 
HCA1 (6) 
HCA1(7) 
HCA2 (1) 
HCA2 (2) 
HCA2 (3) 
HCA2 (4) 
HCA2(5) 
HCA2 (6) 
HCA2 (7) 
HCA3(1) 
HCA3 (2) 
HCA3 (3) 
HCA3 (4) 
HCA3(5) 
HCA3 (6) 


0641 
0277 
.0314 
.0003 
0002 
0002 
.0031 
.1149 
.0499 
0559 
.0040 
0.0007 
. 0007 
. 0046 
. 0632 
. 0266 
0292 
0027 
0007 
0007 


of 


HCA3(7) = -0.0021 
HCA4(1) = 0.0670 
HCA4(2) = -0.0283 
HCA4(3) = 0.0310 
HCA4(4) = -0.0012 
HCA4(5) = -0.0008 
HCA4(6) = 0.0008 
HCA4(7) = 0.0031 
HCA5(1) = 0.0732 
HCA5(2) = -0.0301 
HCA5(3) = 0.0316 
HCA5(4) = -0.0045 
HCA5(5) = -0.0016 
HCA5(6) = 0.0016 
HCA5(7) = 0.0024 
HCA6(1) = -0.0263 
HCA6(2) = -0.0056 
HCA6(3) = -0.0091 
HCA6(4) = 0.0006 
HCA6(5) = -0.0144 
HCA6(6) = 0.0144 
HCA6(7) = -0.0013 
HCA7(1) = -0.5796 
HCA7(2) = -1.6357 
HCA7(3) = -0.0880 
HCA7(4) = -0.1590 
HCA7(5) = -1.8067 
HCA7(6) = 1.8067 
HCA7(7) = -0.0808 


Hydrodynamic coefficient prediction equation 


C1 = 8.023E-03 
RATIO(1) = 0.5686 
RATIO(2) = -1.4357 
RATIO(3) = ~Q.2658 
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RATIO(4) = 0.2675 
RATIO(5) = 1.1781 
RATIO(6) = -30.5114 
RATIO(7) = 0.8149 


DO 5000 I=1,7 
HC (1I)=CHCA1 (1) * FN**2+HCA2 (I) *FN*FM 
& +HCA3 (I) *FM**2+HCA4 (I) *FN 
+HCA5 (I) *FM+HCA6 (1) 
+HCA7 (I) * (VOLUME/ (L*L*L) -C1)) /RATIO(T) 


5000 CONTINUE © 


HC (8) = ~-6.33E-04 

ZQDOT = ~§.33E-04*0.5*RHO*L**4 
HC (8) = ZQDOT 

ZWDOT =  HC(5)*0.5*RHO*L**3 

ZQ = HC (3) *0 .5*RHO*#L**3 
ZW = HC(1)*0.5*RHO*L**2 
MQDOT = HC(7)*0.5*RHO*L**5 

MWDOT = HC(6)*0.5*RHO*L**4 

MQ = HC(4)*0.5*RHO*L**4 
MW = HC(2)*0.5*RHO*L**3 
RATIO1 = 0.92943 

LY = IY/RATIOI 
WRITE(*,4001) IY 

NDi = 0.5*RHO*L**2 

ZGB = ZG-ZB 


DEFINE THE LENGTH X AND BREADTH B TERMS FOR THE INTEGRATION 


DO 333 I=0,21 
XL(I+1)= I*LB/21.0 
BR(I+1) =DIAM*XL(I+1)/LB 
333 CONTINUE 
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DO 334 I=1,2 
XL(22+I)= LB+I*LM/2.0 
BR(22+1)=DIAM 
334 CONTINUE 
DO 335 I=1,30 
E WRITE(*,*) I 
XL(1I+24)= XLCI+23)+1./4.*(L-XL(1I+23)) 
IF ((C(XL(1I+24)-LB-LM) **2/(LN**2)) .GT.1.0) THEN 
BR(I+24)=0.0 
ELSE 
BR(I+24) =DIAM*SQRT (1.0-( (XLC1+24) ~LB-LM) **2/ (LN**2) ) ) 
ENDIF 
335 CONTINUE 
XL(55) 
BR(55) 0 
DO 102 N = 1,55 
XLC(N) = XLCN)-L+LCB 
102 CONTINUE 





WRITE(20,7001) XL 
WRITE(20,7001) BR 


DO 104 K = 1,55 
VECO (K) =BR (K) 
VEC1 (K) =XL(K) *BR(K) 
VEC2 (K) =XL (K) *XL (K) *BR(K) 
VEC3 (K) =XL (K) *XL (K) #XL (K) *BR(K) 
VEC4 (K) =XL (K) *XL(K) *XL (K) *XL (K) *BR (CK) 
104 CONTINUE 
CALL TRAP(55,VECO,XL,EO0) 
CALL TRAP(55,VEC1,XL,E1) 
CALL TRAP(55,VEC2,XL,E2) 
CALL TRAP(55,VEC3,XL,E3) 
CALL TRAP(55,VEC4,XL,E4) 


EPSILON = 0.001 
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XGMIN=-0.01 
XGMAX=+0.0O1 
IXG=80 
XGMIN=XGMIN*L 
XGMAX=XGMAX*L 


C(aaemaemcssoesossansecsursessesesseesesceseesensessesesssssasesae=s 
DO 1 IT=1,1XG 
C WRITE (*,3001) IT,IXG 
XG=XGMIN+ (XGMAX-XGMIN) * (IT-1) / (IXG-1) 
XGB=XG-XB 
DV=(MASS-ZWDOT) * (IY-MQDOT) 
& ~ (MASS*XG+ZQDOT) * (MASS*XG+MWDOT) 


CD6=CD/ (3.DO*EPSILON*DV) 

DW1=CD6* ( (IY-MQDOT) * (-E0) + (MASS*XG+ZQD0T) *E1) 
DW2=CD6* ( (IY-MQDOT) * (3*E1) - (MASS*XG+ZQDOT) *3*E2) 
DW3=CD6* ( (IY-MQDOT) * (~3*E2) +(MASS*XG+ZQDOT) *3+*E3) 
DW4=CD6* ( (IY-MQDOT) * (E3) - (MASS*XG+ZQDOT) *E4) 

DQ1=CD6* ( (MASS-ZWDOT) * (E1) +(MASS*XG+MWDOT) * (-E0O) ) 
DQ2=CD6* ( (MASS-ZWDOT) * (-3*E2) +(MASS*XG+MWDOT) * (3+E1)) 
DQ3=CD6* ( (MASS-ZWDOT) * (3*E3) + (MASS*XG+MWDOT) * (-3*E2) ) 
DQ4=CD6* ( (MASS-ZWDOT) * (-E4) + (MASS*XG+MWDOT) * (E3) ) 
THETAO=ATAN (-XGB/ZGB) 

AAO=(MASS-ZWDOT) * (IY-MQDOT) 


& — (MWDOT+MASS*XG) * (ZQDOT+MASS*XG) 
BBO= (-ZWDOT*MASS-MASS*MW-ZQ*MASS ) *XG 

& +(-MASS*MQ+ZWDOT*MQ-ZQDOT*MW 

& -ZQ*MWDOT-MASS*MWDOT-IY*ZW+MQDOT*ZW) 


CCO=-MASS*ZW*XG+MQ* ZW-ZQ*MW-MASS*MW 
CCi=(-MASS*XG+ZWDOT*XG+MASS*XB-ZWDOT*XB) *SIN (THETAO) 


& + (-MASS*ZB-ZWDOT*ZG+ZWDOT*ZB+MASS*ZG) *COS (THETAO) 
DD1= (ZW*XG-ZW*XB) *SIN (THETAO) + (ZW*ZB-ZW*ZG) 
& *COS (THETAO) 


C After applying AD=BC ( Routh Criterion ), we manage to calculate 
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C the critical speed UO. 


WEI=BBO*CCO/ (AAO*DD1-BBO*CC1 ) 
UO=DSQRT (1556 .2363/WEI) 

U=U0 

WRITE (*,*) U/8.0,XG/L 


DETERMINE [A] AND [B] COEFFICIENTS 


A11DV=(CIY-MQDOT) *ZW+ (MASS*XG+ZQDOT) *MW 
A12DV=(IY-MQDOT) * (MASS+ZQ) + (MASS*XG+ZQDOT) * (MQ-MASS*XG) 
A13DV=- (MASS*XG+ZQDOT) *WEIGHT 
A21DV=(MASS-ZWDOT) *MW+ (MASS*XG+MWDOT) *ZW 
A22DV=(MASS-ZWDOT) * (MQ-MASS*XG) + (MASS*XG+MWDOT) * (MASS+ZQ) 
A23DV=- (MASS-ZWDOT ) *WEIGHT 


A11=A11DV/DV 
A12=A12DV/DV 
A13=A13DV/DV 
A21=A21DV/DV 
A22=A22DV/DV 
A23=A23DV/DV 


C11DV=CIY-MQDOT) *MASS*ZG 
C12DV=- (MASS*XG+ZQDOT ) *MASS*ZG 
C21DV=- (MASS-ZWDOT) *MASS*ZG 
C22DV= (MASS*XG+MWDOT) *MASS*ZG 


C1i1=C1iiDV/DV 
C1i2=C1i2DV/DV 
C21=C21DV/DV 
C22=C22DV/DV 


eee ee fen a Dc Me et eins SD EE SS A i a A A i a Sr A A a cin i a ee Se miki eee SD ewe fee See See See me melee see eee 
eS SS Se A SY A SS Ged i Sy ee ee ee Se ee ee ee ee ee ee oe 


EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 


K1i=- (XGB*SIN (THETAO) -ZGB*COS (THETAO) ) 
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12 
11 


14 


13 


16 


K2=-(1./6.)*(ZGB*COS (THETAO) -XGB*SIN (THETAO) ) 


A(1,1)=0.0 
A(1,2)=0.0 
A(1,3)=1.0 
A(2,1)=A13*K1 
A(2,2)=A11*U 
A(2,3)=A12*U 
A(3,1)=A23*K1 
A(3,2)=A21*U 
A(3,3)=A22*U 
DO 11 I=1,3 
DO 12 J=1,3 
ASAVE(I, J)=ACI, J) 
CONTINUE 
CONTINUE 


CALL RG(3,3,A,WR,WI,1,YYY,IV1,FVi, IERR) 
CALL DSOMEGCIEV,WR,WI,OMEGA, CHECK) 


WRITE (*,*) IEV 


WRITE (*,*) (WRCIWR) , IWR=1,3) 
WRITE (*,*) (WICIWI) ,IWI=1,3) 


OMEGAO=OMEGA 
DO 5 I=1,3 
T(I,1)= YYY(I, IEV) 
T(I,2)=-YYY(1I, IEV+1) 
CONTINUE 
IF (IEV.EQ.1) GO TO 13 
IF (IEV.EQ.2) GO TO 14 
STOP 3004 
DO 6 I[=1,3 
T(I,3)=YYY(I, 1) 
CONTINUE 
GO TO 17 
DO 16 I=1,3 
T(I,3)=YYY(I, 3) 
CONTINUE 
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CONTINUE 


NORMALIZATION OF THE CRITICAL EIGENVECTOR 


INORM=1 
IF (INORM.NE.O) CALL NORMAL(T) 


INVERT TRANSFORMATION MATRIX 


DO 2 I=1,3 
DO 3 J=1,3 
TINV(I,J)=0.0 
TSAVEC(I, J)=TC(I, J) 
CONTINUE 
CONTINUE 
CALL DLUD(3,3,TSAVE,3,TLUD,IVLUD) 
DO 4 I=1,3 
IF (IVLUD(I).EQ.0) STOP 3003 
CONTINUE 
CALL DILU(3,3,TLUD,IVLUD,SVLUD) 
DO 8 I=1,3 
DO 9 J=1,3 
TINV(I, J)=TLUD(I, J) 
CONTINUE 
CONTINUE 


CHECK Inv(T) *A*T 


IMULT=1 

IF (IMULT.EQ.1) CALL MULTCTINV, ASAVE,T,A2) 
IF (IMULT.EQ.0) STOP 

P=A2(3,3) 

PEIG=P 

WRITE (*,4001) (A2(1,JA2) ,JA2=1,3) 

WRITE (*,4001) (A2(2,JA2) ,JA2=1,3) 

WRITE (*,4001) (A2(3,JA2) ,JA2=1,3) 
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PAUSE 


DEFINITION OF Ni; 


N11=TINV (1,1) 
N21=TINV(2,1) 
N31=TINV (3, 1) 
Ni2=TINV(1,2) 
N22=TINV (2,2) 
N32=TINV (3,2) 
Ni3=TINV(1,3) 
N23=TINV (2,3) 
N33=TINV (3,3) 


DEFINITION OF Mij - 


M11=T(1,1) 
M21=T(2,1) 
M31=T(3, 1) 
M12=T(1,2) 
M22=T (2,2) 
M32=T(3, 2) 
M13=T(1,3) 
M23=T(2,3) 
M33=T (3,3) 


DEFINITION OF Lij 


L25=C11*M31*M31+C12*M21*M31 
L26=2*C11*M31*M32+C12* (M21+*M32+M22+*M31) 
L27=C11*M32*M32+C12*M22*M33 
L35=C22*M31*M31+C21*M21*M31 
L36=2*C22*M31*M32+C21* (M21*M32+M22*M31) 
L37=C22*M32*M32+C21*M33*M22 


DEFINITION OF ALFA, BETA, GAMA 
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C 
C 
C 
C 
& 
C 
& 
& 
C 
& 
& 
C 
C 
C 


Di =N32*L25 + N33*L35 
D2 =N32*L26 + N33*L36 
D3 =N32*L27 + N33*L37 


Dii=-P 
D12=0MEGAO 
D21=-2*0OMEGAO 
D22=-P 
D23=2*0MEGAO 
D32=-OMEGAO 
D33=-P 


BETA=(D2-D21*D1/D11-D23*D3/D33) 

/ (D22-D21*D12/D11-D23*D32/D33) 
ALFA=(D1i-D12*BETA) /D11 
GAMA=(D3-D32*BETA) /D33 


L21A=2*C11*ALFA*M31*M33+C12*ALFA 
* (M21*M334+M23*M31) 


L22A=2*C11*ALFA*M32*M33 + 2*C11*BETA*M31*M33 
+ C12*ALFA* (M22*M33+M32*M23) 
+ (C12*BETA* (M21*M33+M23*M31) 
L23A=2*C11*GAMA*M31*M33+2*C11*BETA*M32*M33 
+ C12*GAMA* (M21*M33+M23*M31) 


+ C1i2*BETA* (M22*M33+M23*M32) 


L24A=2*C11*GAMA*M32*M33+C12*GAMA 
* (M22*M334+M23*M32) 


L31A=2*C22*ALFA*M31*M33+C21*ALFA 
* (M21#*M334+M23*M31) 


L32A=2*C22*ALFA*M32*M33+2*C22*BETA*M31*M33 
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+ C21*ALFA*(M22*M33+M32*M23) 
+ C21*BETA*(M21+*M33+M23*M31) 


L33A=2*C22*GAMA*M31*M33+2*C22*BETA*M32*M33 
+ (C21*GAMA* (M21*M334+M23*M31) 
+ C21*BETA* (M22*M33+M23*M32) 


L34A=2*C22*GAMA*M32*M33+C21*GAMA 
* (M22*M33+M23*M32) 


L21=L21A+A13*K2*M11**3+DW1*M21**3 
+DW2*M314*M21**2 
+DW3*M21*M31**2+DW4*M31*«3 

L22=L22A+3*A13*K2*M12*M11**2+3*DW1*M22*M21**2 

+DW2* (2*M21*M22*M31+M324*M21**2) 
+DW3* (24M214*M31*M32+M22*M31**2) 
+3*DW4*M32*M31**2 

L23=L23A+3*A13*K2*M11*M12**2+3*DW1*M21*M22**2 

+DW2* (M31*M22**2+2*M21*M22*M32) 


+ DW3* (M21*M32**2+2*M22*M31*M32) 
+ 3*DW44*M31*M32**2 
L24=L24A+A13*K2*M12**3+DW1*M22**3 
+DW2*M32*M22**2 


+DW3*M22*M32**2+DW44M32**3 
L31=L31A+A23*K24*M11**3+DQ14*M21**3 
+DQ2*M31*M21**2 
+DQ3*M21*M31**2+DQ4*M31**3 
L32=L32A+3*A23*K2*M12*M11**2+3*DQ14*M22*M21**2 
+DQ2* (24M21*M22+*M31+M324*M21**2) 
+DQ3%* (2*M21*M31*M32+M22*M31**2) 


+3*DQ4*M32*M31**2 
L33=L33A+3*A23*K2*M11*M12**2+3%*DQ1*M21*M22**«2 
+ DQ2* (M31*M22**2+2*M21*M22*M32) 
+ DQ3* (M21*M32**2+24M22*M314*M32) 
+ 3*DQ4*M31*M32**2 


L34=L34A+A23*K24*M12**3+DQ1*M22**3 
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+DQ2*M32*M22**2 
+DQ3*M22*M32**2+DQ4*M32**3 


R11=N12*L21+N13*L31 
R1i2=N12*L22+N13*L32 
R13=N12*L23+N13*L33 
R14=N12*L24+N13*L34 
R21=N22*L21+N23*L31 
R22=N22*L22+N23*L32 
R23=N22*L23+N23*L33 
R24=N22*L24+N23*L34 


EVALUATE DALPHA AND DOMEGA 


UINC=0 .001 
UR =U+UINC 
UL =U-UINC 
U =UR 


A(1,1)=0.0 
A(1,2)=0.0 
A(1,3)=1.0 
A(2,1)=A13*K1 
A(2,2)=A11*U 
A(2,3)=A12*U 
A(C3,1)=A23*K1i 
AC3,2)=A21*U 
A(3,3)=A22*U 


CALL RG(3,3,A,WR,W1I,0,YYY,IV1,FV1, IERR) 
CALL DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEOS 

OMEGR=FREQ 


--- e 


U=UL 
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A(1,1)=0.0 
A(1,2)=0.0 
A(1,3)=1.0 
A(2,1)=A13*Ki 
A(2,2)=A11*U 
A(2,3)=A12*U 
A(3,1)=A23*K1 
A(3,2)=A21*U 
A(3,3)=A22*U 


CALL RG(3,3,A,WR,WI,0,YYY,IV1,FV1, IERR) 
CALL DSTABL(DEOS, WR,WI, FREQ) 
ALPHL=DEO0S 

OMEGL=FREQ 


DALPHA=(ALPHR-ALPHL) / (UR-UL) 
DOMEGA= (OMEGR-OMEGL) / (UR-UL) 


EVALUATION OF HOPF BIFURCATION COEFFICIENTS 


COEF 1=3 .0*R11+R13+R22+3 .0*R24 
COEF 2=3 .0*R21+R23-R12-3 .0*R14 
AMU2 =-COEF1i/(8.0*DALPHA) 
BETA2=0 .25*COEF1 
TAU2 =- (COEF2-DOMEGA*COEF1/DALPHA) / (8 .0*O0MEGAO) 
PER =2.0*3.1415927/0MEGAO 
PER =PER*U/L 
WRITE (20,2001) XCB 
WRITE (20,2001) XG/L,COEF1 
1 CONTINUE 
8889 CONTINUE 
8888 CONTINUE 
8887 CONTINUE 
8886 CONTINUE 
STOP 
1001 FORMAT (’ ENTER NUMBER OF DATA LINES’) 


69 


1002 FORMAT (’ ENTER UO, ZG, AND DSAT’) 
1003 FORMAT (’ ENTER BOW PLANE TO STERN PLANE RATIO’) 
1004 FORMAT (’ ENTER ZG’) 
2001 FORMAT (2E14.5) 
4001 FORMAT (1F15.5) 
7001 FORMAT (6F15.5) 
3001 FORMAT (215) 
END 
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